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We present an algorithm for finding the probabilities of rare events in nonequilibrium processes. 
The algorithm consists of evolving the system with a modified dynamics for which the required 
event occurs more frequently. By keeping track of the relative weight of phase-space trajectories 
generated by the modified and the original dynamics one can obtain the required probabilities. The 
algorithm is tested on two model systems of steady-state particle and heat transport where we find 
a huge improvement from direct simulation methods. 
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I. INTRODUCTION 

A rare event is one which occurs with a very small 
probability. However, when they do occur they can have 
a huge effect and so it is often important to estimate the 
actual probability of their occurrence. Examples where 
rare events are important are in banking and insurance, 
in biological systems where important processes such as 
genetic switching and mutations occur with extremely 
small rates, and in nucleation processes. Rare events are 
also of importance in nonequilibrium processes such as 
charge and heat transport in small devices and transport 
in biological cells. The functioning of nano-electronic de- 
vices can be affected by rare large-current fluctuations 
and it is important to know how often they occur. 

In this paper our interest is in predicting probabilities 
of rare fluctuations in transport processes. A number of 
interesting results have been obtained recently on large 
fluctuations away from typical behavior in nonequilib- 
rium systems. These include results such as the fluctua- 
tion theorems [3-0| and the Jarzynski relation Q . In the 
context of transport one typically considers an observ- 
able, say Q, such as the total number of particles or heat 
transferred across an object with an applied chemical po- 
tential or temperature difference. This is a stochastic 
variable and for a given observation time r this will have 
a distribution P(Q,r). The various general results that 
have been obtained for P(Q, r) give some quantitative 
measure of the probability of rare fluctuations. Analytic 
computations of the tails of P(Q,t) for any system are 
usually difficult. This is also true in experiments or in 
computer simulations since the generation of rare events 
requires a large number of trials. 

For large r the probabilities of large fluctuations show 
scaling behavior P(Q,t) ~ e - T /(Q/ T ) i where the func- 
tion f(q) is known as the large deviation function [jjllioj. 
For a few model systems exact results have been ob- 
tained 043 for either f(q) or its Legendre transform 
/i(A), which can be defined in terms of the character- 
istic function as /i(A) = lim T _ ! . 00 T -1 ln(e _A( ^). Recently 
an algorithm has been proposed to compute /i(A). 
However, as has been pointed out in Ref. 12] there may 



be problems in obtaining the tails of /i(A) using the al- 
gorithm of Ref. [Til] - The algorithm proposed in this 
paper is complementary to the one discussed in Ref. [llj 
in the sense that we obtain P(Q,t) directly. Our algo- 
rithm, based on the idea of importance sampling, com- 
putes P(Q,t) for any given r and accurately reproduces 
the tails of the distribution. Algorithms based on impor- 
tance sampling |13[ have earlier been used in the study 
of equilibrium systems and in the study of tran- 



sition rate processes |16l - ll8j . However, we are not aware 



of any applications to the study of large fluctuations of 
currents in nonequilibrium systems and this is the main 
focus of this paper. Here we choose two prototype models 
of transport, namely, heat conduction across a harmonic 
chain and particle transport in the symmetric simple ex- 
clusion process. We illustrate the implementation of im- 
portance sampling in the computation of large fluctua- 
tions of currents in these two nonequilibrium systems. 

Consider a system with a time evolution described 
by the stochastic process x(t). For simplicity we as- 
sume for now that x(t) is an integer- valued variable and 
time is discrete. Let us denote a particular path in 
configuration space over a time period r by the vector 
x(t) := {x(t)\t — 1, 2, . . . , r} and let Q be an observable 
which is a function of the path x(t). We will be inter- 
ested in finding the probability distribution P(Q, r) of Q 
and especially in computing the probability of large devi- 
ations about the mean value (Q) . As a simple illustrative 
example consider the tossing of a fair coin. For r = N 
tosses we have a discrete stochastic process described by 
the time series x(N) — {xi] where Xi = 1 if the out- 
come in the ith trial is heads and Xi = — 1 otherwise. 
Suppose we want to find the probability of generating 
Q heads (thus Q = YliLi ^x 4 ,i)- An example of a rare 
event is, for example, the event Q — N. The probability 
of this is 2~ N and if we were to simulate the coin toss 
experiment we would need more than 2 N repeats of the 
experiment to realize this event with sufficient frequency 
to calculate the probability reliably. For large N this 
is clearly very difficult. The importance sampling algo- 
rithm is useful in such situations. The basic idea is to 
increase the occurrence of the rare events by introducing 
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a bias in the dynamics. The rare events are produced 
with a new probability corresponding to the bias. How- 
ever, by keeping track of the relative weights of trajec- 
tories of the unbiased and biased processes it is possible 
to recover the required probability corresponding to the 
required unbiased process. 

II. THE ALGORITHM 

We now describe the algorithm in the context of evalu- 
ating P(Q, t) for the stochastic process x(r). We denote 
the probability of a particular trajectory by 'P(x) . By 
definition: 

P(Q,t)=^ q>q(x) P(x). (1) 

X 

For the same system let us consider a biased dynamics 
for which the probability of the same path x is given by 
■Pb(x). Then we have: 

P(g,r)^^^ Q(x) e^( x )n(x), (2) 

X 

Thus in terms of the biased dynamics, P(Q, r) is the av- 
erage (SQ^Q( x )e~ w )b and in a simulation we estimate this 
by performing averages over M realizations to obtain: 

P e (Q,T) = ±J25 QMxr) e- w ^) , (4) 

r 

where x r denotes the path for the rth realization. For 
M — > oo we obtain P e (Q,r) — > P(Q,t) which is the re- 
quired probability. Note that the weight factor W is a 
function of the path. In a simulation we know the details 
of the microscopic dynamics for both the biased and un- 
biased processes. Thus we can evaluate W for every path 
x generated by the biased dynamics. A necessary require- 
ment of the biased dynamics is that the distribution of 
Q that it produces [i.e., Pb(Q,r) — (<5Q ! Q( x ))b] should be 
peaked around the desired values of Q for which we want 
an accurate measurement of P(Q,r). As we will see the 
required dynamics can often be guessed from physical 
considerations. 

We first explain the algorithm for the coin tossing ex- 
periment. In this case we consider a biased dynamics 
where the probability of heads is p and that of tails 
is 1 — p. If we take p « 1 then the event Q = N, 
which was earlier rare, is now generated with increased 
frequency and we can use Eq. Q to estimate the re- 
quired probability P(Q = N,N). For any path con- 
sisting of Q heads the weight factor is simply given by 
e~ w = {l/2) N /\p Q {l-p) N - Q ]. Choosing p = 0.95 it 
is easy to see that for N = 100 we can get the required 
probability P(Q = N,N) with more than 1% accuracy 
using only M — 10 7 realizations as opposed to at least 



M = 10 30 required by the unbiased dynamics. Note that 
for this example W has the same value for all paths with 
the same Q. In general of course W depends on the de- 
tails of the path, e.g. for a random walk with a waiting 
probability. We will now illustrate the algorithm with 
non-trivial examples of computing large deviations in two 
well known models in nonequilibrium physics. These are 
the (i) symmetric simple exclusion process (SSEP) with 
open boundaries and (ii) heat conduction across a har- 
monic system connected to Langevin reservoirs. 



III. SYMMETRIC SIMPLE EXCLUSION 
PROCESS 



This is a well studied example of an interacting 
stochastic system consisting of particles diffusing on a 
lattice with the constraint that each site can have at most 
one particle. Here we restrict ourselves to one-dimension 
and study the case of an open system where a linear chain 
with L sites is connected to particle reservoirs at the two 
ends. The dynamics can be specified by the following 
rules: (a) a particle at any site I = 1,2, ... ,L can jump 
to a neighboring empty site with unit rate; (b) at I = 1 a 
particle can enter the system with rate a (if it is empty) 
and leave with rate 7. At site N a particle can leave or 
enter the system with rates /3 and 5, respectively. The 
biased dynamics can be realized in various ways, for ex- 
ample, by introducing asymmetry in the bulk hopping 
rates or by changing the boundary hopping rates. 

For SSEP, the configuration of the system at any time 
is specified by the set C = n\{t), r^i(t)} where 

ni(t) (0 or 1) gives the occupancy of the Ith. site. The dy- 
namical rules specify the matrix element W(C,C) giving 
the transition rate from configuration C to C. We write 
W{C,C) = Wi + W-i + Wo where Wi and W_i corre- 
spond to transitions whereby a particle enters the system 
from the left bath or leaves the system into the left bath, 
respectively, while Wo corresponds to all other transi- 
tions. At long times the system will reach a steady state 
with particles flowing across the system and we are here 
interested in the current fluctuations in the wire. Specifi- 
cally, let Q be the net particle transfer from the left reser- 
voir into the system during a time interval r. For a fixed 
r we want to obtain the distribution P(Q, r) of Q, in the 
steady state of the system. It is useful to define the joint 
probability distribution function R{Q,C, t) for Q number 
of particles transported and for the system to be in state 
C, given that at r = the system is in the steady state. 
Clearly P(Q,t) = J2c R (Q> C > T )- We also define the 
characteristic functions R(z,C,t) = P(Q, C, t)z q 

and P(z, r) = R(z,C,t). It is then easy to obtain 
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FIG. 1. (Color online) Plot of P(Q) for r = 15 for the one- 
site SSEP model with a = /3 = 3.0, 7 = 5 = 3.0. MC refers 
to direct Monte Carlo simulations. Left bias corresponds to 
a' = ft = 3.8,7' = 6' = 2.2 and right bias to a' = ft = 
2.2,7' = 5' = 3.8. 



FIG. 2. (Color online) Plot of P(Q) for r = 15 for the three- 
site SSEP model with a = /9 = 4.0, 7 = 5 = 2.0. MC refers 
to direct Monte Carlo simulations. For left (right) bias simu- 
lations, the particles in bulk hop to the left (right) with rate 
4 and to the right (left) with unit rate. The boundary rates 
are kept unchanged. 



the following master equation 



dR{z, C, t) 



Y,[zWi{C,C')+W (C,C) 



+ ~w_i(e,c')J r(z,c',t) 



(5) 



The general solution of this equation for arbitrary L is 
difficult but for L = 1 an explicit solution can be ob- 
tained for R(z,C',t) and P(z, r). We will here first 
discuss a special case a — /3 — 7 = S for which 
P(z, r) can be inverted explicitly. The choice of steady 
state initial conditions gives the solution: P(Q, r) = 
(e- 2Q 72)[/ 2 Q_ 1 (2ar) + 2I 2Q (2a T ) + J 2Q+1 (2ar)]. In 
Fig. ([T]) we plot the exact distribution along with a di- 
rect simulation of the above process with averaging over 
5 x 10 8 realizations. As we can see the direct simulation 
is accurate only for events with probabilities of O(10 -8 ). 
Now we illustrate our algorithm using a biased dynamics. 
We consider biasing obtained by changing the boundary 
transition rates. We denote the rates of the biased dy- 
namics by a', (3', 7', 8' and these are chosen such that 
Pb(Q) bas a peak in the required region. In our simula- 
tion we consider a discrete-time implementation of SSEP. 
For every realization of the process over a time r (after 
throwing away transients) the weight factor W is dynam- 
ically evaluated. For instance, every time a particle hops 
into the system from the left reservoir, W is incremented 
by — In (a /a'). In Fig. ([I]) we see the result of using our 
algorithm with two different biases. Using the same num- 
ber of realizations we are now able to find probabilities 
up to O(10~ 16 ) and the comparison with the exact result 
is excellent. 

We next study the case with L — 3 with rates chosen 
such that the system reaches a nonequilibrium steady 
state with (Q) > 0. Finding R(z,C,t) analytically in- 



volves diagonalizing an 8 x 8 matrix. We do this numeri- 
cally and after an inverse Laplace transform find P(Q, t). 
In Fig. ([2]) we show the numerical and direct simulation 
results for this case and also the results obtained using 
the biased dynamics; in this case we consider a biased 
dynamics with asymmetric bulk hopping rates. Again 
we find that the biasing algorithm significantly improves 
the accuracy of finding probabilities of rare events using 
the same number of realizations (5 x 10 8 ). 



IV. HEAT CONDUCTION 

Next we consider the problem of heat conduction 
across a system connected to heat reservoirs modeled by 
Langevin white-noise reservoirs. Here we are interested 
in the distribution of the net heat transfer Q from the 
left bath into the system over time t. First let us con- 
sider the simple example of a single Brownian particle 
connected to two baths at temperatures Ti and T 2 . This 
model was studied recently by Visco @ who obtained an 
exact expression for the characteristic function of Q. The 
equation of motion for the system is given by: 



-(71 + 72)« + V2A m + V2D 



2 V2 



(6) 



where 771,2 are Gaussian delta-correlated noises with zero 
mean and unit variance , thus (?)i{t)T)j(t')) = 6ijS(t — t') 
and Di — TiTj. The heat flow from the left bath into 
the system in time r is given by Q(t) = Jiv 2 + 

\[TD\r\\v) dt. For the single Brownian particle in this 
problem it is sufficient to specify the state by the veloc- 
ity v(t) alone. If we choose T\ > T 2 then P(Q, r) will 
have a peak at Q > 0. It is clear that to use the bi- 
asing algorithm to compute probabilities of rare events 
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FIG. 3. (Color online) Plot of P(Q) for r = 200 for heat 
conduction across a single free particle with 71 = 0.8, 72 = 
0.2, Tj = 1.1875, T 2 = 0.25. The parameters have been cho- 
sen to correspond to a region in parameter space where the 
fluctuation theorem is not satisfied Q. MC refers to di- 
rect Monte Carlo simulations. The left bias corresponds to 
71 = 71.72 = 72/20, T{ = 21, n = 20T 2 . 

with Q < we can choose a biased dynamics with tem- 
peratures of left and right reservoirs taken to be T{ and 
T2 with T[ < T' 2 . The calculation of the weight fac- 
tor W is somewhat tricky since computing P [«(£)] from 
V[i]i(t), 372 (i)] is non-trivial. Also one cannot eliminate 
771 to express Q as a functional of only the path v. To 
get around this problem we note the following mapping 
of the single-particle system to the over-damped dynam- 
ics of two coupled oscillators given by the equa- 
tions of motion: x\ = —71(^1 — X2) + \/1D\T]\ , £2 = 
—72(22 — %i) — V2£ ) 2?72. The variable x\ — X2 = £12 
satisfies the same equation as v in Eq. Thus with 
the same definition for Q as given earlier we can use the 
above equations for x\ and X2 to find P(Q,t). In this 
case we do not have the problem as earlier and both Q 
and W can be readily expressed in terms of {x\, 22}. Let 
us denote by 7-, T(, D\ the parameters of the biased sys- 
tem. Also let r][ 2 be the noise realizations in the biased 
process that result in the same path {xi, X2} as produced 
by ?7:l.2 for the original process. Choosing Di = D\ for 
i = 1, 2 it can be shown that: 

W= f dt [(7??/2 + % 2 /2)-(r ? f/2 + ^/2)]. (7) 
Jo 

Using the equations of motion we can express 771,2, ?7i 2 
in terms of the phase-space variables and this gives: 

1 f T 

+ J dt[2(j 2 ~ 12)^2x12 + {ll - 72)^12]) 

Q = I dtx\X\2- 
Jo 
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FIG. 4. (Color online) Plot of P(Q) for r = 100 for heat con- 
duction across two particles connected by a harmonic spring 
with unit spring constant and 71 = 72 = \/2, 21 = 10, T 2 = 
12. MC refers to direct Monte Carlo simulations. The left 
bias corresponds to j[ =71,72 = 72/2, T{ = T 1: Ti = 2T 2 
and right bias to j[ = 7 i/2, j' 2 = 72, T{ = 2Ti , T 2 = T 2 . 

Thus W and Q are easily evaluated in the simulation 
using the biased dynamics. In Fig. ((3]) we show results 
for P(Q,t) obtained both directly and using the biased 
dynamics. Again we see that for the same number of 
realizations (10 9 ) one can obtain probabilities about 10 8 
times smaller than using direct simulations. The compar- 
ison with the numerical results obtained from the exact 
expression for (e~ A( 2) Q also shows the accuracy of the 
algorithm. 

It is easy to apply the algorithm to more complicated 
cases. For example consider a one-dimensional chain of 
L particles connected to heat reservoirs at the two ends 
with the following equations of motion: 

mivi = fi + <J/,i[-7iUi + \J2D\ 771] 

+ 6i, L [-72V L + V2IhV2}, l = l,2,...,JV, (8) 

where // = —d Xl U and U({xi}) is the potential energy of 
the system. The net heat transfer from the left bath into 
the system is given by Q = jl (—71^1 + \/2D 1T71U1) and 
using Eqs. ((5J) this can be expressed in terms of {xi,v{\ 
as Q = J Q r dtv\ {m\V\ — To apply our algorithm we 
consider a biased dynamics where the Hamiltonian evo- 
lution is unchanged while the bath dynamics has new 
parameters 7^ , 72 , T-[ , which are chosen so that Pb{Q) 
has a peak in the required region. Choosing D' t — Di we 
again find W by using Eqs. ((SJ) in Eq. ([7]), as for the single 
particle case. Thus both Q and W can be expressed in 
terms of the path and so are readily evaluated for every 
realization of the biased dynamics. 

As an example we study the case L = 2 with U = 
(xi — X2) 2 /2 and with mi = ni2 = 1. For the special 
parameters 71 = 72 = \[2 we use the results in Ref. Q 
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to obtain (e- AC ?) - e^ T with //(A) = \/2{l - [1 + 

~ ^)] 1/6 } • This can be inverted to numer- 
ically compute P(Q, t) at large r. In Fig. (Q| we give the 
comparison between the analytical distribution and that 
obtained by the biasing method. 

V. CONCLUSION 

In conclusion, we have presented an algorithm for com- 
puting the probabilities of rare events in various nonequi- 
librium processes. The algorithm is an application of 
importance sampling and consists in using a biased dy- 
namics to generate the required rare events. This al- 
gorithm is straightforward to understand and also to 
implement. The error in the estimate of P(Q, t) is 
« (e- 2W 6 Q , Q S b /2 /[MPbiQ)} 1 / 2 . In the systems that we 



have studied we find that the error can be made small by 
choosing the biased dynamics carefully. We have applied 
the algorithm to two different models of particle and heat 
transport and shown that in both cases it gives excellent 
results. We note, however, that, in general, the fluctua- 
tions in W grow with r and with the system size, hence 
the errors are large and finding an appropriate biased 
dynamics is not always easy. Further work is necessary 
for improving the efficiency of the algorithm for general 
systems. 
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